--- title: patrones_espacio_temporales keywords: fastai sidebar: home_sidebar summary: "Métodos para explorar los patrones espacio-temporales de la delincuencia y su evolución" description: "Métodos para explorar los patrones espacio-temporales de la delincuencia y su evolución" nb_path: "01_hotspots.ipynb" ---
{% raw %}
{% endraw %} {% raw %}
{% endraw %}

kde2D

construye_malla

{% raw %}

construye_malla[source]

construye_malla(datos, size)

Regresa una malla (np.meshgrid) ajustada al extent de los datos, con el tamaño de celda especificado.

Args: datos (GeoDataFrame): carpetas o víctimas size (float): tamaño de las celdas (en las unidades de la proyección)

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
carpetas = get_carpetas_from_api(10000)
carpetas = carpetas.to_crs(32614)
malla = construye_malla(carpetas, 200)
{% endraw %}

ajusta_bandwidth_kde

{% raw %}

ajusta_bandwidth_kde[source]

ajusta_bandwidth_kde(datos, bandwidth_space, size=1000, malla=None, n_jobs=-1, metric='euclidean')

Regresa el valor de bandwidth con mejor log likelihood.

Parametros:

datos (GeoDataFrame):  víctimas o carpetas
bandwith_space (np.linspace):  con el espacio de búsqueda
size (float): Tamaño de la celda (en las unidades de la proyección).
              Si se especifica malla se ignora
malla (np.meshgrid): la malla en la que se va a ajustar el KDE, si es None se calcula
n_jobs (int): número de procesos a usar (default = -1)
metric (str): métrica a usar para calcular las distancias (default euclidean)
{% endraw %} {% raw %} 
def ajusta_bandwidth_kde(datos, bandwidth_space, size=1000, 
                         malla=None, n_jobs=-1, metric="euclidean"):
    """ Regresa el valor de bandwidth con mejor log likelihood.

        Parametros:
        
            datos (GeoDataFrame):  víctimas o carpetas
            bandwith_space (np.linspace):  con el espacio de búsqueda
            size (float): Tamaño de la celda (en las unidades de la proyección).
                          Si se especifica malla se ignora
            malla (np.meshgrid): la malla en la que se va a ajustar el KDE, si es None se calcula
            n_jobs (int): número de procesos a usar (default = -1)
            metric (str): métrica a usar para calcular las distancias (default euclidean)
    """
    if malla is None:
        xx, yy = construye_malla(datos, size)
    else:
        xx = malla[0]
        yy = malla[1]        
    xy_sample = np.vstack([yy.ravel(), xx.ravel()]).T
    x = datos.geometry.x.to_numpy()
    y = datos.geometry.y.to_numpy()
    xy_train  = np.vstack([y, x]).T
    grid = GridSearchCV(KernelDensity(metric=metric), bandwidth_space, n_jobs=n_jobs)
    grid.fit(xy_train)
    return grid.best_estimator_.bandwidth
{% endraw %}

Esta función ajusta el ancho de banda para calcular un KDE en 2 dimensiones. Puede ser muy tardado.

El ancho de banda se puede calcular en cualquyier sistema de coordenadas, sin embargo es conveniente proyectarlos a coordenadas planas.

Primero ajustamos especificando el tamaño de la malla

{% raw %} 
params = {'bandwidth': np.linspace(10, 10000, 100)}
bw = ajusta_bandwidth_kde(carpetas, params, size=1000)
print(bw)

312.72727272727275
{% endraw %}

También podemos especificar la malla primero

{% raw %} 
params = {'bandwidth': np.linspace(10, 10000, 100)}
bw = ajusta_bandwidth_kde(carpetas, params, malla=malla)
print(bw)

312.72727272727275
{% endraw %}

kde2D

Regresa la superficie del KDE con los parámetros especificados

{% raw %}

kde2D[source]

kde2D(datos, bandwidth, size=1000, malla=None)

Regresa una matriz con la densidad de kernel para los datos.

Parametros:

datos (GeoDataFrame):  víctimas o carpetas
bandwith: ancho del kernel gaussiano
size (float): Tamaño de la celda (en las unidades de la proyección).
              Si se especifica malla se ignora
metric (str): métrica a usar para calcular las distancias (default euclidean)
malla (np.meshgrid): la malla en la que se va a ajustar el KDE, si es None se calcula
{% endraw %} {% raw %}
{% endraw %}

Si no mandamos una malla como argumento, se calcula utilizando size

{% raw %} 
xx, yy, zz = kde2D(carpetas, bw, size=1000)
fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection='3d')
ax = ax.plot_surface(xx, yy, zz,cmap='viridis', edgecolor='none')
{% endraw %}

También se puede especificar una malla con el tamaño deseado y calcular el KDE

{% raw %} 
malla = construye_malla(carpetas, 1000)
xx, yy, zz = kde2D(carpetas, bw, malla=malla)
fig = plt.figure(figsize=(10,10))
ax = plt.axes(projection='3d')
ax = ax.plot_surface(xx, yy, zz,cmap='viridis', edgecolor='none')
{% endraw %}

malla_comun

{% raw %}

malla_comun[source]

malla_comun(datos, size)

Regresa las mallas (formato np y xarray) más pequeñas (de tamaño size) que contienen a todos los datos.

Args: datos(list(GeoDataGrames)): lista con los datos para la malla común. size (float): tamaño de las celdas en la malla.

{% endraw %} {% raw %}
{% endraw %}

Esta función se usa dentro de serie_tiempo_kde_categoria y no se llama directamente.

get_lista_datos

{% raw %}

get_lista_datos[source]

get_lista_datos(carpetas, fechas, categorias, offset)

Regresa una lista de GeoDataFrames con los datos segmentados en fechas para la categoría seleccionada.

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
carpetas = get_historico_carpetas()
carpetas = agregar_categorias_carpetas(carpetas)
carpetas = carpetas.to_crs(32614)
fechas = pd.date_range(start='1/1/2020', end='3/1/2020', freq='M').to_list()
categorias = ["Homicidios dolosos"]
datos = get_lista_datos(carpetas, fechas, categorias, "30 days")
datos[0].head()
ao_hechos mes_hechos fecha_hechos ao_inicio mes_inicio fecha_inicio delito fiscalia agencia unidad_investigacion ... calle_hechos calle_hechos2 colonia_hechos alcaldia_hechos competencia longitud latitud tempo geometry categoria
895799 2020.0 Enero 2020-01-01 04:50:00 2020 Enero 2020-01-01 05:59:10 HOMICIDIO POR ARMA DE FUEGO FISCALÍA DE INVESTIGACIÓN ESTRATÉGICA DEL DELI... 2 1 CON DETENIDO 1 C/D ... 15 Y AVENIDA GUADALUPE NaN GUADALUPE PROLETARIA GUSTAVO A MADERO NaN -99.156937 19.521753 NaN POINT (483534.561 2158567.439) Homicidios dolosos
895810 2020.0 Enero 2020-01-01 07:30:00 2020 Enero 2020-01-01 08:43:13 HOMICIDIO POR ARMA DE FUEGO FISCALÍA DE INVESTIGACIÓN ESTRATÉGICA DEL DELI... 2 1 CON DETENIDO 1 C/D ... PIRULES NaN TLALPEXCO GUSTAVO A MADERO NaN -99.128763 19.544280 NaN POINT (486492.381 2161057.661) Homicidios dolosos
895835 2020.0 Enero 2020-01-01 11:00:00 2020 Enero 2020-01-01 11:04:10 HOMICIDIO POR ARMA DE FUEGO FISCALÍA DE INVESTIGACIÓN TERRITORIAL EN IZTAP... IZP-9 UI-3SD ... elisa NaN 2A. AMPLIACIÓN SANTIAGO ACAHUALTEPEC IZTAPALAPA NaN -99.015140 19.347795 NaN POINT (498409.849 2139310.995) Homicidios dolosos
895836 2020.0 Enero 2020-01-01 09:55:00 2020 Enero 2020-01-01 11:05:21 HOMICIDIO POR ARMA BLANCA FISCALÍA DE INVESTIGACIÓN TERRITORIAL EN MIGUE... H1 UI-1SD ... F.C. NACIONALES DE MEXICO AND. NUEVA GALICIA NUEVA ESPAÑA AZCAPOTZALCO NaN -99.189655 19.498750 NaN POINT (480099.007 2156025.577) Homicidios dolosos
896408 2020.0 Enero 2020-01-02 19:20:00 2020 Enero 2020-01-02 19:53:43 HOMICIDIO POR ARMA DE FUEGO FISCALÍA DE INVESTIGACIÓN ESTRATÉGICA DEL DELI... 2 1 CON DETENIDO 1 C/D ... 25 DE SEPTIEMBRE DE 1873 PUENTE QUEMADO LEYES DE REFORMA 3A SECCIÓN IZTAPALAPA NaN -99.061941 19.376333 NaN POINT (493495.517 2142469.932) Homicidios dolosos

5 rows × 21 columns

{% endraw %}

KDE_hotspots

Regresa la serie de tiempo de KDEs para las categorías seleccionadas. La serie de tiempo se construye a partir de los datos de incidentes (víctimas o carpetas de investigación) y se agrega en intervalos de tiempo definidos por el parámetro fechas. En este parámetro deben venir los extremos derechos de los intervalos (regulares) para agregar y debe coincidir con el parámetro offset que se usa internamente para calcular el primer intervalo de agregación.

Cuando esta función se usa para calcular diréctamente los Hotspots pára una categoría, no es necesario llamarla con una malla, la función la va a calcular a partir de los datos, lo que es necesario es fijar el parámetro size que indica el tamaño de los pixeles en la malla resultante.

Para usar la función es necesario que labase de incidentes contenga la columna categoria.

{% raw %}

KDE_hotspots[source]

KDE_hotspots(carpetas, fechas, size, categorias, offset, malla=None, bw=1000)

Ajusta kdes egregando los datos sobre cada categoria e intervalo de fecha.

Args: carpetas (GeoDataFrame): Incidentes (carpetas/victimas), si viene malla se ignora size size (float): Tamaño de la celda (en las unidades de la proyección) fechas list(pd.datetime): extremos (derechos) de los intervalos de tiempo categorias: Lista de categorías para calcular el KDE offset: intervalo para agregar antes de la primera fecha, p.ej: "30 days" si los intervalos son mensuales malla (dict): Obtenida con malla_comun. bw: bandwidth para la estimación del KDE

Salida: serie_kde: xr.DataArray con coordenadas (lat, long, tiempo) que contiene los arrays para cda intervalo.

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
fechas = pd.date_range(start='1/1/2019', end='1/1/2020', freq='M').to_list()
categorias = ["Homicidios dolosos"]
xr_kde = KDE_hotspots(carpetas, fechas, 1000, categorias, "30 days", bw=1000)
assert type(xr_kde) == xr.DataArray
xr_kde

<xarray.DataArray 'Densidad' (tiempo: 12, latitud: 49, longitud: 40)>
array([[[2.35342140e-95, 3.26405373e-91, 3.45287329e-86, ...,
         1.88889308e-31, 7.81712655e-31, 1.19606204e-30],
        [4.00050431e-88, 5.54122870e-84, 2.82477511e-80, ...,
         2.20508720e-27, 9.06597278e-27, 1.38057662e-26],
        [2.50170065e-81, 3.46518689e-77, 1.76572679e-73, ...,
         9.55127104e-24, 3.89229718e-23, 5.88894101e-23],
        ...,
        [1.46696807e-53, 6.25631743e-48, 9.81920184e-43, ...,
         2.08730206e-37, 6.20890950e-42, 6.86910197e-47],
        [9.03945700e-56, 3.85738497e-50, 6.05559721e-45, ...,
         5.59847363e-40, 1.65382275e-44, 1.79985312e-49],
        [2.05108754e-58, 8.75281404e-53, 1.37409616e-47, ...,
         5.54831898e-43, 1.63760305e-47, 1.77849903e-52]],

       [[1.31792911e-67, 5.16965706e-62, 7.45996748e-57, ...,
         1.15844838e-20, 5.83565178e-22, 1.08388389e-23],
        [9.71536379e-64, 3.81091050e-58, 5.49925617e-53, ...,
         5.80492128e-18, 2.92501581e-19, 5.43297480e-21],
        [2.63470261e-60, 1.03347811e-54, 1.49133938e-49, ...,
         1.07016109e-15, 5.39357035e-17, 1.00184070e-18],
...
        [1.10833679e-57, 1.53175687e-52, 7.78788308e-48, ...,
         3.59329966e-43, 3.43703619e-48, 1.25577806e-53],
        [1.57812688e-61, 2.18103896e-56, 1.10896499e-51, ...,
         7.52353944e-46, 6.49518009e-51, 2.21073550e-56],
        [8.26645978e-66, 1.14250136e-60, 5.81065084e-56, ...,
         7.23322001e-49, 5.41923247e-54, 1.63332504e-59]],

       [[1.42046716e-86, 2.12065796e-82, 1.16471479e-78, ...,
         3.52322398e-22, 1.58040023e-24, 2.60795104e-27],
        [6.88541162e-80, 1.02790194e-75, 5.64529891e-72, ...,
         7.41900817e-20, 3.32791849e-22, 5.49167756e-25],
        [1.22803789e-73, 1.83315221e-69, 1.00671514e-65, ...,
         5.74720864e-18, 2.57800523e-20, 4.25418277e-23],
        ...,
        [3.97735036e-59, 6.56257079e-54, 3.98345278e-49, ...,
         7.84263318e-39, 7.08556625e-43, 2.35501454e-47],
        [5.02419020e-63, 8.28984147e-58, 5.03189912e-53, ...,
         2.44791433e-42, 2.21159548e-46, 7.35062752e-51],
        [2.33476844e-67, 3.85233444e-62, 2.33835283e-57, ...,
         2.81091709e-46, 2.53947236e-50, 8.44035803e-55]]])
Coordinates:
  * latitud   (latitud) float64 2.116e+06 2.117e+06 ... 2.163e+06 2.164e+06
  * longitud  (longitud) float64 4.65e+05 4.66e+05 ... 5.03e+05 5.04e+05
  * tiempo    (tiempo) datetime64[ns] 2019-01-31 2019-02-28 ... 2019-12-31
{% endraw %}

La función regresa un xr.DataArray de xarray con los KDEs etiquetados por los extremos de los intervalos. Podemos utilizar geoviews para geerar fácilmente una visualización de los resultados.

{% raw %} 
%%output  size=200 widget_location='bottom'
xr_dataset = gv.Dataset(series, 
                        kdims=["tiempo", "longitud", "latitud"],                          
                        crs=crs.UTM(zone=14))
(xr_dataset
 .to(gv.Image, ["longitud", "latitud"])
 .opts(alpha=0.7) * gv.tile_sources.CartoDark())
{% endraw %}

En el ejemplo anterior seleccionamos por el índice (xr_kde.isel(tiempo=0)), también podemos seleccionar por el valor de la etiqueta

razones_de_incidentes

Esta función calcula la setrie de tiempo de la razón de una categoría de incidentes con respecto a todas las demás. Se usa internamente para calcular más adelante los mapas de intensidad y significancia.serie_tiempo_kde_categoria

{% raw %}

razones_de_incidentes[source]

razones_de_incidentes(carpetas, fechas, categorias, offset, size, bw)

Regresa un xr.Dataset con los xr.DataArray para lass series por categoría, base y las razones.

{% endraw %} {% raw %}
{% endraw %}

Calculamos la serie de mapas de razones

{% raw %} 
fechas = pd.date_range(start='1/1/2019', end='1/1/2020', freq='M').to_list()
categorias = ["Homicidios dolosos"]
series = razones_de_incidentes(carpetas,fechas, ["Homicidios dolosos"], "30 days", 1000, 1000)
series

<xarray.Dataset>
Dimensions:           (latitud: 51, longitud: 42, tiempo: 12)
Coordinates:
  * latitud           (latitud) float64 2.115e+06 2.116e+06 ... 2.165e+06
  * longitud          (longitud) float64 4.642e+05 4.652e+05 ... 5.052e+05
  * tiempo            (tiempo) datetime64[ns] 2019-01-31 ... 2019-12-31
Data variables:
    Serie categorias  (tiempo, latitud, longitud) float64 5.225e-105 ... 4.23...
    Serie base        (tiempo, latitud, longitud) float64 4.968e-72 ... 1.161...
    Serie razones     (tiempo, latitud, longitud) float64 1.052e-33 ... 3.65e-05
{% endraw %}

El resultado es un xr.Dataset que contiene las dos series (base y de categoría). Con geoviews podemos, por ejemplo, comparar el resultado de KDE_hotspots para la categoría especificada con la serie para todas las demás categorías.

{% raw %} 
%%output  size=200 widget_location='bottom'
xr_dataset = gv.Dataset(series, 
                        kdims=["tiempo", "longitud", "latitud"],
                        vdims=["Serie categorias"],
                        crs=crs.UTM(zone=14))
(xr_dataset
 .to(gv.Image, ["longitud", "latitud"], ["Serie categorias"])
 .opts(alpha=0.7) * gv.tile_sources.CartoDark())
{% endraw %}

Podemos visualizar las tres superficies resultantes

smer_KDE

Esta función calcula los Hotspots basados en la idea de SMER (Social Media Events Ratio). El resultado es un xr.Dataset con:

  • Los Hotspots KDE para la categoría seleccionada
  • Los Hotspots KDE para el agregado de todas las demás categorías
  • La serie de tiempo de las razones entre los hotspots para la categoría seleccionada contra las demás
  • La intensidad (z-scores) de la serie de razones
  • Una medida de la significancia estadística de la intensidad

La idea es que esta función nos permite ver qué tanto resalta la categoría seleccionada con respecto al resto de las categorías y cómo se va moviendo en el tiempo. Es una especie de medida del saliency de la categoría de interés y nos permite identificar lugares en donde esta categŕía empieza a desplazar a las demás.

{% raw %}

smer_KDE[source]

smer_KDE(carpetas, fechas, categorias, offset, size, bw)

Regresa los mapas de razon y las intensidades de la categoría para las fechas seleccionadas.

{% endraw %} {% raw %}
{% endraw %}

La función regresa un xr.Dataset con variables para cada una de las series.

{% raw %} 
fechas = pd.date_range(start='1/1/2019', end='1/1/2020', freq='M').to_list()
categorias = ["Homicidios dolosos"]
series = smer_KDE(carpetas, fechas, ["Homicidios dolosos"], "30 days", 1000, 1000)
series

<xarray.Dataset>
Dimensions:           (latitud: 51, longitud: 42, tiempo: 12)
Coordinates:
  * latitud           (latitud) float64 2.115e+06 2.116e+06 ... 2.165e+06
  * longitud          (longitud) float64 4.642e+05 4.652e+05 ... 5.052e+05
  * tiempo            (tiempo) datetime64[ns] 2019-01-31 ... 2019-12-31
Data variables:
    Serie categorias  (tiempo, latitud, longitud) float64 5.225e-105 ... 4.23...
    Serie base        (tiempo, latitud, longitud) float64 4.968e-72 ... 1.161...
    Serie razones     (tiempo, latitud, longitud) float64 1.052e-33 ... 3.65e-05
    intensidad        (tiempo, latitud, longitud) float64 -0.3018 ... -0.3054
    p_values          (tiempo, latitud, longitud) float64 0.8462 ... 0.5385
{% endraw %}

En general, las variables importantes de esta función son la Intensidad y la significancia de los hotspots identificados con SMER. Podemos ver mapas de ambos de forma fácil usando geoviews

{% raw %} 
%%output widget_location='bottom', size=150
xr_dataset = gv.Dataset(series, 
                        kdims=["tiempo", "longitud", "latitud"],
                        vdims=["intensidad", "p_values"],
                        crs=crs.UTM(zone=14))
int_gv = (xr_dataset
         .to(gv.Image, ["longitud", "latitud"], ["intensidad"])
         .opts(alpha=0.7) * gv.tile_sources.CartoDark())
p_gv = (xr_dataset
        .to(gv.Image, ["longitud", "latitud"], ["p_values"])
        .opts(alpha=0.7) * gv.tile_sources.CartoDark())
lo = int_gv + p_gv
lo
{% endraw %}